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Abstract: The current trend of the civil aviation technology is to modernize the legacy air 
traffic control (ATC) system that is mainly supported by many ground based navigation 
aids to be the new air traffic management (ATM) system that is enabled by global 
positioning system (GPS) technology. Due to the low receiving power of GPS signal, it is a 
major concern to aviation authorities that the operation of the ATM system might 
experience service interruption when the GPS signal is jammed by either intentional or 
unintentional radio-frequency interference. To maintain the normal operation of the ATM 
system during the period of GPS outage, the use of the current radar system is proposed in 
this paper. However, the tracking performance of the current radar system could not meet 
the required performance of the ATM system, and an enhanced tracking algorithm, the 
interacting multiple model and probabilistic data association filter (IMMPDAF), is 
therefore developed to support the navigation and surveillance services of the ATM 
system. The conventional radar tracking algorithm, the nearest neighbor Kalman filter 
(NNKF), is used as the baseline to evaluate the proposed radar tracking algorithm, and the 
real flight data is used to validate the IMMPDAF algorithm. As shown in the results, the 
proposed IMMPDAF algorithm could enhance the tracking performance of the current 
aviation radar system and meets the required performance of the new ATM system. Thus, 
the current radar system with the IMMPDAF algorithm could be used as an alternative 
system to continue aviation navigation and surveillance services of the ATM system during 
GPS outage periods. 
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1. Introduction 

According to the International Civil Aviation Organization's plan, the current technology trend is to 
develop and implement Communications, Navigation, Surveillance, and Air Traffic Management 
(CNS/ATM) systems based on the Global Positioning System (GPS) to replace legacy Air Traffic 
Control (ATC) systems based on ground-based radar [1]. The Next Generation Air Transportation 
System (NextGen) will modernize the air traffic control system in the United States, and many of the 
foundational elements required to meet the predicted capacity and efficiency improvements rely on 
widespread use of precision Positioning, Navigation, and Timing (PNT) services provided by the 
GPS [2]. The Federal Aviation Administration (FAA) needs to ensure a sufficient backup PNT 
capability to reduce the risks of aviation users if GPS becomes unavailable. As described in the FAA 
Flight Plan [2], the air traffic system of the future will be much more dependent on GPS. This strategic 
goal will not be achievable without a NextGen Alternate PNT (APNT), especially in the event of a 
GPS loss of service. The APNT program ensures alternate PNT services provided by the ATC system, 
minimizes the economic impact of GPS outages, and supports air transportation's timing needs. The 
existing legacy navigation and surveillance infrastructure, i.e., VHF Omnidirectional radio Range 
(VOR), Distance Measuring Equipment (DME), Non-Directional Beacon (NDB), and secondary radar, 
may not achieve the desired level of performance for NextGen operation [3]. 

To improve ATC tracking performance, one solution is to upgrade the ATC radar beacon system. 
However, this solution is expensive. A more cost-effective solution is to improve the ATC tracking 
algorithms, which does not require the implementation of new radar facilities. The present study thus 
proposes a radar tracking algorithm that meets the requirements of APNT systems. The major challenge 
and difficulty for ATC tracking is tracking target maneuvering in a cluttered environment [4-7]. 
Bar-Shalom and Blom [4] introduce a tracking algorithm called the Interacting Multiple Model (IMM) 
estimator, which provides tracking estimates with significant noise reduction and fast response to 
sequences of aircraft maneuver modes [4,5]. However, the tracking of an aircraft in a cluttered 
environment might be a challenge due to the several observations for a single aircraft under such 
environment [6,7]. That is, some tracking measurements do not originate from the target aircraft. 
Therefore, the present study utilizes the Probabilistic Data Association (PDA) filter [6,7] to assign 
weights to the validated measurements. The PDA filter can extend the tracking capability to a highly 
cluttered environment. In order to gain possible improvement on the tracking performance, this study 
combines the IMM estimator and PDA filters to create an IMMPDA filter (IMMPDAF). Most studies 
on radar tracking algorithms use simulated data for analysis. To further evaluate the IMMPDAF, this 
study uses real flight radar data collected by the Civil Aeronautics Administration (CAA) of Taiwan. A 
filter that uses the nearest neighbor method and the standard Kalman filter (NNKF) is commonly used 
for radar tracking systems. Therefore, this paper compares the tracking performance of the IMMPDAF 
and the NNKF in terms of positioning accuracy. Finally, the validation of the IMMPDAF under APNT 
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requirements is discussed. Due to the complexity of the IMMPDAF, this paper also conducts a 
computational load study on the IMMPDAF and the NNKF. 

The rest of this paper is organized as follows: Section 2 introduces the algorithms of the IMM 
estimator and the PDA filter. Then, the IMM estimator and the PDA filter are integrated. The dynamic 
model setting of the IMMPDAF for the ATC scenario is also discussed in Section 2. The 
implementation of the IMMPDAF to a radar system is discussed in Section 3. The tracking 
performance of the IMMPDAF is compared with that of the NNKF. In addition, the computation loads 
of these filters are evaluated. Section 4 concludes this work. 



2. IMMPDAF 

2.1. IMM Estimator 



The IMM estimator is a good candidate for a radar tracking system [4]. The IMM estimator 
provides significant noise reduction and a fast response. Previous investigations have found that the 
IMM estimator is a cost-effective technique for tracking a maneuvering target [4]. In the IMM 
estimator, several Kalman filters are used in parallel. Each Kalman filter uses a different dynamic 
model. The IMM estimator can be separated into four steps. In the first step, the state estimate and 
covariance are mixed. These mixed estimates and covariances are the inputs of the Kalman filter. The 
equations are: 

u lJ (k-^) = ^ LA — - m 
x° j (k-i)=±x;(k-i)u iJ (k-i) (2) 

i=l 

^(*-l) = ±^(*-l){^(*-l) + [^(*-l)-^(*-l)][^(*-l)-^(*-0] r ) (3) 
1=1 L J 

In Equations (1-3), the subscripts and / indicate different Kalman filters. There are up to r 
different dynamic models of Kalman filter, the total number for m„ X] , X] , P° and P* is r, and the total 
number for Uy and py is r 2 . w, (k—\) represents the model probabilities in the previous time 
(k—l). Since the initial model probabilities have little impacts on the results over time, they can be 
conveniently chosen. In this study, the initial model probabilities are set to be rectangular functions, py 
is an element of the Markov transition matrix; there is probability py that the target will transit from 
model i to model j at each time step. Uy represents the conditional probability of the target in the j 
model state, which transited from the i model state. According to [5], the IMM algorithm performance 
is robust to the choice of transition matrix. Xt(k~l) and P + t (k—\) are the updated state and covariance, 
respectively, from the Kalman filter in the previous time step (k~\). X] (k-\) and Pj (k~\) are the 
mixed state and the mixed covariance, respectively, for each Kalman filter as inputs at time k. Then, 
the second step implements the multiple models of the Kalman filter: 
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j X-{k) = F,X«{k-l) 
\P7(k) = F j .P°(k-l).F; + Q j 

K j = P j (A) /// (// /' (*)•//,'+*,)' 

*;(*)=*;(*)+M z (*)-^- x ;(*)) 

i?(*) = (/-*,.^K(*) 



(4) 



(5) 



The details of Kalman filter can be found in [8]. In Equations (4,5) the subscript j indicates it is the j 
th Kalman filter, and Equations (4,5) are applied for each Kalman filter. In Equation (4), Fj is the 
dynamic model and Qj is the covariance of the process noise. XJ (k) and PJ (k) are the uncorrected 
predicted state and covariance, respectively. In Equation (5), Hj is the observation matrix, Rj is the 
covariance of the measurement noise, and Kj is the Kalman gain. Z(k) is the measurement at time k and 
/ is the identity matrix. Xj(k) and Pj(k) are the corrected predicted state and covariance, respectively. 
In the third step, the updated model probabilities are calculated as: 

.V, (*) = // •/> 
y j (k) = Z(k)-H j .X:(k) 

d J 2 (k) = y J (k) T S J (ky l y J (k) 



exp 



M*) = 



-dHk) 



(6) 



Uj (k) 



^S~[k)\ 

1=1 



(7) 



In Equation (6), the subscript j indicates it is the y'th Kalman filter; exp[] means the exponential 
operation and | | represents calculating the determinant of the matrix. The likelihood function A 7 - is the 
probability density function of the estimated measurement Hj -XJ (k), which is a normal distribution 
with mean Z{k) and covariance Sj(k). In Equation (7), the subscripts i and j indicate different Kalman 
filters. Uj(k) represents the latest model probabilities at time k. In the last step, the corrected predicted 
state Xj (k) and P] (k) covariance from the Kalman filters are combined. The weighting factor is the 
latest model probabilities Uj(k): 



x(k) = ^x;(k)- Uj (k) 

./=! 



p(k) = ±u j (k){p;(k) + [x](k)-x(k)lx](k)-x(k)J] 



(8) 



(9) 



X(k) and P(k) are the combinations of the system state and covariance from 1 to rth Kalman filters, 
respectively. Using Equations (1-9), the IMM estimator calculates the system state and covariance 
from time k~ 1 to k. 
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2.2 PDA Filter 

The PDA filter is designed for tracking a target in a cluttered environment based on the Kalman 
filter [6,7]. The PDA filter obtains an estimator which incorporates all the returning measurements that 
might originate from the target of interest rather than select only one of them. The resulting estimator 
sets a gate to determine whether the measurements are valid. The association probabilities of the target 
being tracked for each validated measurement are then calculated. These probabilities are assigned to 
all of the validated measurements with different weights. Next, the validated measurements are 
combined with different weights based on location, and then the combined measurement is used to 
update the state estimate of the target. The PDA filter can extend the tracking capability to a highly 
cluttered environment. 

The probabilities are used for weighting the correction and covariance in the PDA filter. The basic 
assumptions and theories for the PDA filter can be found in [6,7]. Here, a brief introduction of its 
functions and equations is given. Because the PDA filter is based on the Kalman filter, the first 
equations are the same as those for the Kalman filter: 



Similarly, the first step is to predict the system state and covariance. In Equation (10), the subscript j 
indicates it is the y'th PDA filter, X*(k~l) is the calculated state at time k—\, F is the dynamic model, 
and X] (k) is the predicted state at time k . H is the observation matrix and z. (k) is the predicted 

measurement at time k. Pj(k—l) is the calculated covariance at time k—l and Qj is the process noise 
covariance. The second step is to determine which measurements can be used to update the state and 
covariance: 



The main purpose of Equation (11) and (12) is to validate the measurements by Chi-Square test, and 
it includes all the measurements returned by all the radars. In Equation (11), the total number of 
measurements is unknown; the subscript i indicates it is the zth measurement. Z, is the innovation of zth 
measurement for y'th PDA filter, Zj(k) is the zth measurement at time k and R is the covariance of the 
measurement noise. Equation (12) is the validation equation, y is the threshold corresponding to the 
gate probability, y can be obtained from Chi-Square tables for a chosen gate probability. Once the zth 
measurement passes the Chi-Square test in Equation (12), it can be utilized in the rest of the PDA 
filter. Then, the association probabilities of each validated measurement are calculated as: 



X~j{k) = F r x]{k-\) 
z j (k) = H j -X-(k) 

p j {k)=F j .p:{k-\).F T j+Qj 



(10) 



Zi(k) = z i (k)-z J (k) 
S j (k) = H j .p;(k).H T j+ R J 

Z/(k) T Sj(k)-Z/(k)<y 



(11) 



(12) 
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r y (*) = (4*/3).^. > /[s^ 

Tjf*_ N [ Z '( k )M k )'> S j( k )} P D 

i{)= m/Vj(k) ^ 



(13) 



I'M 



i-p p 

1 1 D 1 G 



,i = l,---m 



J = 0 



(14) 



In Equation (13), the subscript j indicates it is the y'th PDA filter, V{k) is the volume of the validation 
region, m is number of total validation measurements, and iV[z ( . (k);z(k);S (k)^ represents the normal 
distribution, centered at z(k) with variance S(k). L\ is the likelihood of the zth measurement for y'th PDA 
filter, p, is the associated probability of the zth measurement for y'th PDA filter. In Equations (13) and 
(14), P D and Pg are the target detection probability and the gate probability, respectively. f} x {k) is the 
association probability of the validated measurements. Note that Po(k) is the probability that all 
measurements are not valid [7]. The suggested values of Pd and Pg can be found in [5]. Pd is required 
to be larger than 90% for primary surveillance radar sensor and 98% for secondary surveillance radar 
sensor. P G is the probability that the true measurement from the target is detected and validated; and 
Pg is usually chosen to be at least 95%. If Pg is set too low, the gate y will be too small that no 
measurement can pass the validation and it could be end up losing track. If P G is set too high on the 
other hand, the gate y will be too large that causes too much false alarms (i.e., not true measurement 
from the target) and reduction on the accuracy [5]. The last step is to update the system state and 
covariance. In this step, only the validated measurements and their association probabilities are used: 



v,(*) = £#-Z/(*) 



i=i 



W j (k) = P(k)-H T r S j (k)- 1 

p;(k)=p;(k)-w j (k)-s j (k)-w j (k) 7 



(15) 



m 

X^-Z/^)-Z/(^-v 7 (^).v 7 (/c) 7 



X + j (k) = X-(k) + W j (k).v j (k) 

P : {k)=pi.p- {k) + (x-pi).p;{k) +P ;{k) 



(16) 



In Equations (15,16), the subscript j indicates it is the y'th PDA filter, v(k) is combined innovation; 
W(k) is Kalman gain; P"(k) and P"\k) are the covariances of state updated. In Equation (16), X*~(k) and 
P + (k) are the calculated system state and covariance, respectively at time k. The propagation of the 
PDA filter is shown in Equations (10-16). 
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2.3. IMMPDA Filter 



The IMM estimator and PDA filter can be combined as shown in Figure 1 . The Kalman filters in 
the IMM estimator are replaced by the PDA filters. At the top, Equations (1-3) are implemented to 
deal with the data from the previous epoch. Then, the PDA filters [Equations (10-16)] use the data 
from Equations (1-3) as the initial values and generate the calculated states and covariances. The next 
step is to use Equations (6,7) to obtain the latest model probabilities. The last step is to utilize 
Equations (8,9) to combine the states and covariances with different weights, which are related to the 
latest model probabilities. The combinations of states and covariances, X(k) and P(k) are the outputs of 
the IMMPDAF. In Figure 1, the IMMPDAF has three different models of PDA filter. The number of 
models depends on the number of dynamic motions of the target. In this paper, the tracking target is a 
commercial airline aircraft. The selection of suitable motion models for an ATC system is discussed in 
the next section. 



Figure 1. Flow chart of the IMMPDAF with three models. 
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2.4. Tuning of IMMPDAF 



As described above, the IMMPDAF can be considered as a modification of the Kalman filter. The 
Kalman filter is tuned by changing the process noise covariance Q matrices and the measurement noise 
covariance R matrices. The matrices Q and R adjust the Kalman gain. The R matrices represent the 
accuracy of the measurement device; larger R matrices mean that the measurements can be trusted less. In 
most practical situations, matrices Q and R are unknowns. In this paper, a run-time estimate method [9] is 
used to obtain the measurement noise covariance R. The adaptive measurement noise covariance is 
estimated in Equation (17). v(n)is the innovation in time epoch n, obtained from Equation (15); and 
v avg (k) is average innovation value from time 1 to time k: 



,avg 



Zv,(«) 

(k)=*±— - 



Rj(k) = ^ 



t^j{n)-v; g {k)\ Vj {n)-v^{k)\ 



(17) 



From the conclusions of [9], the run- time estimate method can handle changes of measurement 
noise covariance, and converge faster than that of a standard Kalman filter. 

2.5. Dynamic Models of Commercial Airline Aircraft 

The key to successful target tracking is selecting suitable dynamic models that fit the target 
motion [10,11]. Ideally, each motion should have a correct dynamic model. In an ATC scenario, the 
motions of commercial airline aircraft can be summarized as ascent, descent, en-route, turning, and 
change of speed. These motions can be described as three basic dynamics: constant velocity motion, 
accelerated motion, and turning motion [11,12]. Following [1 1-15] and our previous work [16], this study 
uses the constant velocity model, the acceleration model, and the turning rate estimator model, 
respectively given in Equations (18-20), respectively, for the IMMPDAF. For ATC radar tracking, the 
measurement is the target 3 -dimensional position (x,y, z) of the tracking target in WGS-84 coordinates. In 
Equations (18-20), X is the system state, x, y, z are the positions, x, y, z are the velocities, and a> is the 
turning rate. F represents the dynamic model and t is the time interval. Generally speaking, the time epoch 
is close to 5 seconds. Q is the covariance of process noise and a is zero-mean Gaussian white noise. The a 
value is tuned by experience and testing, and the value is related to the measurement noise covariance: 

x 

y 

z 

X 

y " ~ (18) 

z 



X = 



,F = 



f tf 
0 /, 



Q = 



0 tL 
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(19) 



(20) 



3. Experimental Results and Discussion 

The experiment data are from the CAA of Taiwan. The data formats include radar data (ASTERIX 
AC048) and ADS-B (ASTERIX AC021). The radar data were recorded from local time 12:00:00 
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midnight on 28 March 2012. Hundreds of commercial airline aircraft were tracked and recorded. Two 
cases are used here to investigate whether the performance of the IMMPDAF significantly varies in the 
area surrounding Taiwan. The first case target was from east head to Taiwan and the second case target 
flight away from Taiwan to north. The radar data was processed using the IMMPDAF and the NNKF. 
The ADS-B data have the target 'Callsign' and GPS calculated positions; these were set as the initial 
values of the two methods and regarded as the true positions. As described in the previous section 
IMMPDAF includes three filters of different dynamic models shown in Equations (18-20). The 
parameters of initial model probabilities and transition matrix are shown in Equation (21). The initial 
values of model probabilities have minor effect on the performance, and the IMMPDAF is very robust 
to the choice of transition matrix. Since the values set for P D is 0.9 and Pq is 0.999, and the degree of 
freedom is 3, the gate threshold obtained from Chi-square table is approximately 16. 

W=[V3 1/3 1/3] 
"0.8 0.1 0.1" 
0.1 0.8 0.1 
0.1 0.1 0.8 



(21) 



The performance requirements for APNT can be found in [17]. The minimum performance 
requirements for APNT are shown in Table 1. The FAA categorized the airspace into several zones [17]. 
En-route is the airspace at flight level (FL) 180 to FL 600 (18,000 to 60,000 feet). Terminal is the 
airspace starting at 500 feet above the ground and extending out to 5 statute miles from the airport; it 
then goes up at a 2-degree angle to FL 180. Under 500 feet, the APNT system needs to support the 
LNAV/non-precision approach. Accuracy (95%) is defined as the errors in the horizontal 95% under 
the accuracy requirement. That is, according to a normal distribution, 95% of error (error mean + 2 
standard deviations) has to be under the accuracy requirement. 

Table 1. Required performance for APNT. 





Navigation 


Surveillance 


Accuracy (95%) 


Separation 


NAC P 


En-route 


10 nmi 


5 nmi 


0.1 nmi 


4 nmi 


2 nmi 


Terminal 


1 nmi 


3 nmi 


0.05 nmi 


LNAV 


0.3 nmi 



3.1. Case 1 



In Figure 2, the target heads to Taiwan Taoyuan International Airport (I ATA: TPE, ICAO: RCTP) 
from east to west. The interval between epochs is 5 s. Note that a lot of surveillance radar systems have 
a scan time of close to 5 s. For example, ASR-12 surveillance radar has a scan time of 4 to 6 seconds. 
The tracking performance is shown in Figure 3. In case 1, the target starting at en-route airspace, so the 
accuracy requirement begins from 2 nmi. As shown in Figure 3, the NNKF experiences the instability at 
the beginning and the IMMPDAF performs more stable. The reason of the difference is the strategy to 
deal with the low quality measurements. While the NNKF just chooses to trust one of the 
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measurements and updates, the IMMPDAF gives all the validated measurements weightings and 
updates with all of them. From Figure 3, the tracking performance of both the IMMPDAF and the 
NNKF meets the accuracy requirement. The IMMPDAF outperforms the NNKF, as shown in Table 2. 
The legends for Figures 3 to 5 are illustrated below. 

Figure 2. Flight path for case 1 (target heading toward Taiwan). 
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Figure 3. Performance of IMMPDAF and NNKF for Case 1. 

Case 1 



2.5 



2 



1.5 



I 

LLI 




50 100 150 200 250 300 350 400 
Epoch (5 s) 



--■IMMPDAF error mean + 2 standard deviations 
IMMPDAF error 

" Y^'VWW** li NNKF error mean + 2 standard deviations 

+ NNKF error 
^Accuracy requirement 



3.2. Case 2 



In Case 2, the target takes off from Taiwan Taoyuan International Airport and flies northeast as 
depicted in Figure 4. The tracking performance is shown in the upper plot of Figure 5, and the bottom 
plot of Figure 5 is the zoom-in plot of the upper one which shows the tracking performance for the 
LNAV/non-precision approach region. Note that the NNKF results initially exceed the accuracy 
requirement of APNT whereas the IMMPDAF results remain under the accuracy requirement. That is, 
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the NNKF method does not meet the accuracy requirement (95%) in this case. As a result, the 
IMMPDAF outperforms the NNKF. 

Figure 4. Flight path for Case 2 (target flying away from Taiwan). 
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Figure 5. Performance of IMMPDAF and NNKF for Case 2. 
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The mean accuracies of the horizontal positioning are summarized in Table 2. The IMMPDAF 
outperforms the NNKF for all tested cases. In comparison to the use of the NNKF, the average 
improvement on the positioning (tracking) performance gained from the proposed IMMPDAF is about 
50%. The overall mean accuracies of the horizontal positioning are also listed in Table 2. 
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Table 2. Overall mean accuracies of the horizontal positioning for IMMPDAF and NNKF. 



Accuracy (95%) 


IMMPDAF 


NNKF 


CASES 1 


0.311 nmi 


0.921 nmi 


CASES 2 


0.397 nmi 


0.582 nmi 


AVERAGE 


0.354 nmi 


0.751 nmi 



3.3. Computation Load 

For all the experiments, the computations were performed in MATLAB (R2008a). A PC with an 
Intel Core 2 Duo E7500 2.93-GHz CPU and 2 GB of RAM was used. The computation performance 
was measured using a MATLAB function. The computation performance results are listed in Table 3. 
The values are total seconds that running a case needed. Each test runs 10 times and to obtain the 
average result. The average result is shown that both algorithms in the single target ATC radar tracking 
scenario are below the radar scan time an epoch (5 seconds). Without coding optimization, the 
IMMPDAF requires about nine times than that of the NNKF. 

Table 3. Computation performance of IMMPDAF and NNKF for Cases 1 to 2. Lower 
values are better. 



Cost Time (s) 


IMMPDAF 


NNKF 


CASE 1 


19.754 s 


2.249 s 


CASE 2 


25.110 s 


2.734 s 


AVERAGE 


22.432 s 


2.491 s 



4. Conclusions 

This paper presented a filter based on the Interacting Multiple Model (IMM) estimator and the 
Probabilistic Data Association (PDA) filter, and this filter is thus called IMMPDAF, which was 
applied to the radar tracking system. The real flight radar data was used to evaluate the tracking 
performance of the IMMPDAF, and the conventional Nearest Neighbor Kalman Filter (NNKF) was 
used as the baseline to show the performance gained from the proposed filter. Experimental results 
show that the tracking performance of the IMMPDAF is better than that of the NNKF. In one test case 
(i.e., Case 2), the IMMPDAF meets the accuracy requirement of the Alternate Positioning, Navigation, 
and Timing (APNT), but the NNKF results exceeded the APNT accuracy requirement for some 
periods of time. Although the computation load of the IMMPDAF is nine times that of the NNKF, the 
IMMPDAF is still suitable for ATC radar since processing time of each epoch is below the radar scan 
time (5 s). The proposed IMMPDAF gained about 50% improvement on tacking performance than that 
of the NNKF. Importantly, the radar tracking with the proposed IMMPDAF met the performance 
requirements of the APNT. That is, the current radar with the proposed IMMPDAF is capable to 
provide the navigation and surveillance services for the Air Traffic Management (ATM) system during 
periods of GPS outages. 
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